%% Aim
% Perform advanced statistical analyses on DDM outputs

%% Setup
clc;
clear;
tic;
root_path = pwd; 
addpath(genpath(root_path));

nrole=2;
nprice=10;
nlottery=10;

role_label={'Buyer','Seller'};
role_color_dark={[235/255 112/255 99/255],[106/255 142/255 208/255]}; 
role_color_median={[238/255 137/255 126/255],[144/255 171/255 220/255]}; 
role_color_light={[243/255 171/255 163/255],[165/255 187/255 227/255]};

bias_label={'Val','Res','Val vs. Res'};
bias_color_median={[195/255 205/255 147/255],[192/255 173/255 203/255],[214/255 171/255 128/255],[0.75 0.75 0.75]}; 

%% Load data
beh_data=readtable([root_path,'\data\behavior\behavior.csv']);
sub_data=readtable([root_path,'\data\subject\subject.csv']);
nsub=height(sub_data);

%%










%% Group posteriors identified from HDDM traces

raw_data=readtable([root_path,'\model\hddm\DDM4a\DDM4a_trace.csv']);

ddm_trace_data=table;

ddm_trace_data.t=raw_data.t;

ddm_trace_data.aBuy=raw_data.a_Intercept;
ddm_trace_data.aDelta=raw_data.a_C_Role__T_2_;
ddm_trace_data.aSell=raw_data.a_Intercept + raw_data.a_C_Role__T_2_;

ddm_trace_data.vInterceptBuy=raw_data.v_Intercept;
ddm_trace_data.vInterceptDelta=raw_data.v_C_Role__T_2_;
ddm_trace_data.vInterceptSell=raw_data.v_Intercept + raw_data.v_C_Role__T_2_;

ddm_trace_data.vPrice=raw_data.v_Price;

ddm_trace_data.vLotteryBuy=raw_data.v_Lottery_C_Role__1_;
ddm_trace_data.vLotterySell=raw_data.v_Lottery_C_Role__2_;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);

ddm_trace_data.zBuy=exp(raw_data.z_Intercept_trans)./(1+exp(raw_data.z_Intercept_trans));
ddm_trace_data.zDelta=raw_data.z_C_Role__T_2_;
ddm_trace_data.zSell=ddm_trace_data.zBuy+ddm_trace_data.zDelta;

writetable(ddm_trace_data,[root_path,'\model\hddm\DDM4a\DDM4a_trace_group.csv']);

%% Group posteriors: Descriptive

ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

% likelihood of difference between buyer and seller
ddm.group.likelihood.val=mean(ddm_trace_data.wLotterySell > ddm_trace_data.wLotteryBuy);
ddm.group.likelihood.res=mean(ddm_trace_data.zSell > ddm_trace_data.zBuy);
ddm.group.likelihood.a=mean(ddm_trace_data.aSell > ddm_trace_data.aBuy);
ddm.group.likelihood.zBuy=mean(ddm_trace_data.zBuy <0.5);
ddm.group.likelihood.zSell=mean(ddm_trace_data.zSell >0.5);
ddm.group.likelihood.vInterceptBuy=mean(ddm_trace_data.vInterceptBuy >0);
ddm.group.likelihood.vInterceptSell=mean(ddm_trace_data.vInterceptSell >0);
ddm.group.likelihood.vInterceptDelta=mean(ddm_trace_data.vInterceptDelta > 0);
ddm.group.likelihood.vallog=mean(ddm_trace_data.wLotterySell > ddm_trace_data.wLotteryBuy);
ddm.group.likelihood.aDelta=mean(ddm_trace_data.aDelta > 0);

disp(ddm.group.likelihood);

%% Subject posteriors identified from ***post_subject.csv

ddm_post_data=readtable([root_path,'\model\hddm\DDM4a\DDM4a_post_subject.csv']);

sub_data=[sub_data,ddm_post_data(:,2:end)];
sub_data.ValBias=sub_data.wLotterySell./sub_data.wLotteryBuy;
sub_data.ResBias=sub_data.zDelta;
sub_data.vInterceptBias=sub_data.vInterceptDelta;
sub_data.aBias=sub_data.aDelta;
sub_data.ValBiasLog=log(sub_data.wLotterySell./sub_data.wLotteryBuy);

writetable(sub_data,[root_path,'\data\subject\subject_ddm.csv']);

min(sub_data.ValBiasLog)
max(sub_data.ValBiasLog)
min(sub_data.ResBias)
max(sub_data.ResBias)

% Categorization of participants
sdata=sub_data(sub_data.ValBiasLog>0 & sub_data.ResBias>0,:);
height(sdata)/height(sub_data);

sdata=sub_data(sub_data.ValBiasLog<0 & sub_data.ResBias>0,:);
height(sdata)/height(sub_data);

sdata=sub_data(sub_data.ValBiasLog>0 & sub_data.ResBias<0,:);
height(sdata)/height(sub_data);

sdata=sub_data(sub_data.ValBiasLog<0 & sub_data.ResBias<0,:);
height(sdata)/height(sub_data);

%% Correlations between biases

sdata=sub_data;

[r,p]=corr([sdata.ValBias, sdata.ResBias, sdata.vInterceptBias, sdata.aBias]);
r(2,2)=0;
r(2,3)=0;
r(3,3)=0;
r(1,:)=[];
r(:,4)=[];

p(2,2)=0;
p(2,3)=0;
p(3,3)=0;
p(1,:)=[];
p(:,4)=[];

    figure('Renderer', 'painters', 'Position', [10 10 450 380]);
        subplot(1,1,1);
        A=r;
        h=imagesc(A);
        cmap_up=[linspace(1,192/256,1000)',linspace(1,0,1000)',linspace(1,0,1000)']; 
        cmap_down=[linspace(0/256,1,1000)',linspace(32/256,1,1000)',linspace(102/256,1,1000)']; 
        cmap=[cmap_down;cmap_up];
        colormap(cmap);
        set(gca,'TickDir','out');
        %xlabel('Price');
        %ylabel('Lottery');
        yticks(1:nprice);
        xticks(1:nlottery);
        yticklabels({'2','4','6','8','10','12','14','16','18','20'});
        xticklabels({'1','2','3','4','5','6','7','8','9','10'});
        colorbar;
        caxis([-1 1]);
        hold on
    hold off
    print(1, '-dtiff', [root_path, '\figure\model\ddm\ddm_post_corr.tif.tif'], '-r200');
    close;     

%% Correlations between biases (log)

sdata=sub_data;

[r,p]=corr([sdata.ValBiasLog, sdata.ResBias, sdata.vInterceptBias, sdata.aBias]);
r(2,2)=0;
r(2,3)=0;
r(3,3)=0;
r(1,:)=[];
r(:,4)=[];

p(2,2)=0;
p(2,3)=0;
p(3,3)=0;
p(1,:)=[];
p(:,4)=[];

    figure('Renderer', 'painters', 'Position', [10 10 450 380]);
        subplot(1,1,1);
        A=r;
        h=imagesc(A);
        cmap_up=[linspace(1,192/256,1000)',linspace(1,0,1000)',linspace(1,0,1000)']; 
        cmap_down=[linspace(0/256,1,1000)',linspace(32/256,1,1000)',linspace(102/256,1,1000)']; 
        cmap=[cmap_down;cmap_up];
        colormap(cmap);
        set(gca,'TickDir','out');
        %xlabel('Price');
        %ylabel('Lottery');
        yticks(1:nprice);
        xticks(1:nlottery);
        yticklabels({'2','4','6','8','10','12','14','16','18','20'});
        xticklabels({'1','2','3','4','5','6','7','8','9','10'});
        colorbar;
        caxis([-1 1]);
        hold on
    hold off
    print(1, '-dtiff', [root_path, '\figure\model\ddm\ddm_post_corr_log.tif.tif'], '-r200');
    close;    
    
%% Plot subject posteriors (two)

sdata=sub_data;

figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
x1=sdata.wLotteryBuy;
y1=sdata.zBuy;
[r1,p1]=corr(x1,y1);
scatter(x1,y1,36,role_color_median{1},'filled');
hold on
x2=sdata.wLotterySell;
y2=sdata.zSell;
[r2,p2]=corr(x2,y2);
scatter(x2,y2,36,role_color_median{2},'filled');
ylim([0.3 0.7]);
yticks(0.3:0.1:0.7);
xlim([0 1.2]);
set(ax,'TickDir','out');
ylabel('Starting Point');
xlabel('Relative Weighting');
box on;
line(xlim,[0.5 0.5],'Color',[0 0 0],'LineStyle','--');
print(1, '-dtiff', [root_path,'\figure\model\ddm\ddm_post.tif'], '-r200');
hold off;
close;

%% histogram z
figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
h1=histogram(ddm_trace_data.zBuy);
h1.Normalization='probability';
h1.BinWidth = 0.005;
h1.FaceColor=role_color_light{1};
h1.EdgeColor=role_color_light{1};
hold on
h2=histogram(ddm_trace_data.zSell);
h2.Normalization='probability';
h2.BinWidth = 0.005;
h2.FaceColor=role_color_light{2};
h2.EdgeColor=role_color_light{2};
hold on
xlim([0.3 0.7]);
%ylim([0 800]);

print(1, '-dtiff', [root_path,'\figure\model\ddm\ddm_post_z_distribution.tif'], '-r200');
hold off;
close;

%% histogram w
figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
h1=histogram(ddm_trace_data.wLotteryBuy);
h1.FaceColor=role_color_light{1};
h1.Normalization='probability';
h1.BinWidth = 0.01;
h1.EdgeColor=role_color_light{1};
hold on
h2=histogram(ddm_trace_data.wLotterySell);
h2.FaceColor=role_color_light{2};
h2.Normalization='probability';
h2.BinWidth = 0.01;
h2.EdgeColor=role_color_light{2};
hold on
xlim([0 1.2]);

print(1, '-dtiff', [root_path,'\figure\model\ddm\ddm_post_w_distribution.tif'], '-r200');
hold off;
close;

%% Plot two DDM biases (log)

sdata=sub_data;

figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
x=sdata.ValBiasLog;
y=sdata.ResBias;
[r,p]=corr(x,y);
scatter(x,y,36,[0.75,0.75,0.75],'filled');
hold on
set(ax,'TickDir','out');
ylabel('Response Bias');
xlabel('Valuation Bias');
box on;

line([-1 1.5],[0 0],'Color',bias_color_median{2},'LineStyle','--');
line([0 0],[-0.1 0.2],'Color',bias_color_median{1},'LineStyle','--');

ylim([-0.1 0.2]);
yticks(-0.1:0.1:0.2);

print(1, '-dtiff', [root_path,'\figure\model\ddm\ddm_bias_log.tif'], '-r200');
hold off;
close;

%% histogram val log
figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
h1=histogram(log(ddm_trace_data.wLotterySell./ddm_trace_data.wLotteryBuy));
h1.Normalization='probability';
h1.BinWidth = 0.1;
h1.FaceColor=bias_color_median{1};
h1.EdgeColor=bias_color_median{1};
xlim([-1 1.5]);
%ylim([0 800]);

print(1, '-dtiff', [root_path,'\figure\model\ddm\ddm_val_distribution_log.tif'], '-r200');
hold off;
close;

%% histogram res
figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
h1=histogram(ddm_trace_data.zSell-ddm_trace_data.zBuy);
h1.Normalization='probability';
h1.BinWidth = 0.005;
h1.FaceColor=bias_color_median{2};
h1.EdgeColor=bias_color_median{2};
xlim([-0.1 0.2]);
%ylim([0 800]);

print(1, '-dtiff', [root_path,'\figure\model\ddm\ddm_res_distribution.tif'], '-r200');
hold off;
close;

%%










%% DDM comparison
%dic=readtable([root_path,'\model\hddm\ddm_dic.csv']);

%% DDM2a: controlling valuation bias

raw_data=readtable([root_path,'\model\hddm\DDM2a\DDM2a_trace.csv']);

ddm_trace_data=table;
ddm_trace_data.t=raw_data.t;
ddm_trace_data.aBuy=raw_data.a_Intercept;
ddm_trace_data.aDelta=raw_data.a_C_Role__T_2_;
ddm_trace_data.aSell=raw_data.a_Intercept + raw_data.a_C_Role__T_2_;
ddm_trace_data.vInterceptBuy=raw_data.v_Intercept;
ddm_trace_data.vInterceptDelta=raw_data.v_C_Role__T_2_;
ddm_trace_data.vInterceptSell=raw_data.v_Intercept + raw_data.v_C_Role__T_2_;
ddm_trace_data.vPrice=raw_data.v_Price;
ddm_trace_data.vLotteryBuy=raw_data.v_Lottery;
ddm_trace_data.vLotterySell=raw_data.v_Lottery;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);
ddm_trace_data.zBuy=exp(raw_data.z_Intercept_trans)./(1+exp(raw_data.z_Intercept_trans));
ddm_trace_data.zDelta=raw_data.z_C_Role__T_2_;
ddm_trace_data.zSell=ddm_trace_data.zBuy+ddm_trace_data.zDelta;
ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

%% DDM3a: Controlling response bias

raw_data=readtable([root_path,'\model\hddm\DDM3a\DDM3a_trace.csv']);

ddm_trace_data=table;
ddm_trace_data.t=raw_data.t;
ddm_trace_data.aBuy=raw_data.a_Intercept;
ddm_trace_data.aDelta=raw_data.a_C_Role__T_2_;
ddm_trace_data.aSell=raw_data.a_Intercept + raw_data.a_C_Role__T_2_;
ddm_trace_data.vInterceptBuy=raw_data.v_Intercept;
ddm_trace_data.vInterceptDelta=raw_data.v_C_Role__T_2_;
ddm_trace_data.vInterceptSell=raw_data.v_Intercept + raw_data.v_C_Role__T_2_;
ddm_trace_data.vPrice=raw_data.v_Price;
ddm_trace_data.vLotteryBuy=raw_data.v_Lottery_C_Role__1_;
ddm_trace_data.vLotterySell=raw_data.v_Lottery_C_Role__2_;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);
ddm_trace_data.zBuy=exp(raw_data.z_trans)./(1+exp(raw_data.z_trans));
ddm_trace_data.zDelta(:)=NaN; %raw_data.z_C_Role__T_2_;
ddm_trace_data.zSell=ddm_trace_data.zBuy;%ddm_trace_data.zBuy+ddm_trace_data.zDelta;
ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

%% DDM4: controlling intercept bias

raw_data=readtable([root_path,'\model\hddm\DDM4\DDM4_trace.csv']);

ddm_trace_data=table;
ddm_trace_data.t=raw_data.t;
ddm_trace_data.aBuy=raw_data.a_Intercept;
ddm_trace_data.aDelta=raw_data.a_C_Role__T_2_;
ddm_trace_data.aSell=raw_data.a_Intercept + raw_data.a_C_Role__T_2_;
ddm_trace_data.vInterceptBuy=raw_data.v_Intercept;
ddm_trace_data.vInterceptDelta(:)=NaN;
ddm_trace_data.vInterceptSell=raw_data.v_Intercept;
ddm_trace_data.vPrice=raw_data.v_Price;
ddm_trace_data.vLotteryBuy=raw_data.v_Lottery_C_Role__1_;
ddm_trace_data.vLotterySell=raw_data.v_Lottery_C_Role__2_;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);
ddm_trace_data.zBuy=exp(raw_data.z_Intercept_trans)./(1+exp(raw_data.z_Intercept_trans));
ddm_trace_data.zDelta=raw_data.z_C_Role__T_2_;
ddm_trace_data.zSell=ddm_trace_data.zBuy+ddm_trace_data.zDelta;
ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

%% DDM0: controlling a bias

raw_data=readtable([root_path,'\model\hddm\DDM0a\DDM0a_trace.csv']);

ddm_trace_data=table;
ddm_trace_data.t=raw_data.t;
ddm_trace_data.aBuy=raw_data.a;
ddm_trace_data.aDelta(:)=NaN;
ddm_trace_data.aSell=raw_data.a;
ddm_trace_data.vInterceptBuy=raw_data.v_Intercept;
ddm_trace_data.vInterceptDelta=raw_data.v_C_Role__T_2_;
ddm_trace_data.vInterceptSell=raw_data.v_Intercept + raw_data.v_C_Role__T_2_;
ddm_trace_data.vPrice=raw_data.v_Price;
ddm_trace_data.vLotteryBuy=raw_data.v_Lottery_C_Role__1_;
ddm_trace_data.vLotterySell=raw_data.v_Lottery_C_Role__2_;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);
ddm_trace_data.zBuy=exp(raw_data.z_Intercept_trans)./(1+exp(raw_data.z_Intercept_trans));
ddm_trace_data.zDelta=raw_data.z_C_Role__T_2_;
ddm_trace_data.zSell=ddm_trace_data.zBuy+ddm_trace_data.zDelta;
ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

%% DDM4b: without intercept

raw_data=readtable([root_path,'\model\hddm\DDM4b\DDM4b_trace.csv']);

ddm_trace_data=table;
ddm_trace_data.t=raw_data.t;
ddm_trace_data.aBuy=raw_data.a_Intercept;
ddm_trace_data.aDelta=raw_data.a_C_Role__T_2_;
ddm_trace_data.aSell=raw_data.a_Intercept + raw_data.a_C_Role__T_2_;
ddm_trace_data.vInterceptBuy(:)=0;
ddm_trace_data.vInterceptDelta(:)=0;
ddm_trace_data.vInterceptSell(:)=0;
ddm_trace_data.vPrice=raw_data.v_Price;
ddm_trace_data.vLotteryBuy=raw_data.v_Lottery_C_Role__1_;
ddm_trace_data.vLotterySell=raw_data.v_Lottery_C_Role__2_;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);
ddm_trace_data.zBuy=exp(raw_data.z_Intercept_trans)./(1+exp(raw_data.z_Intercept_trans));
ddm_trace_data.zDelta=raw_data.z_C_Role__T_2_;
ddm_trace_data.zSell=ddm_trace_data.zBuy+ddm_trace_data.zDelta;
ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

%% DDM4ad: without weighting

raw_data=readtable([root_path,'\model\hddm\DDM4ad\DDM4ad_trace.csv']);

ddm_trace_data=table;
ddm_trace_data.t=raw_data.t;
ddm_trace_data.aBuy=raw_data.a_Intercept;
ddm_trace_data.aDelta=raw_data.a_C_Role__T_2_;
ddm_trace_data.aSell=raw_data.a_Intercept + raw_data.a_C_Role__T_2_;
ddm_trace_data.vInterceptBuy=raw_data.v_Intercept;
ddm_trace_data.vInterceptDelta=raw_data.v_C_Role__T_2_;
ddm_trace_data.vInterceptSell=raw_data.v_Intercept + raw_data.v_C_Role__T_2_;
ddm_trace_data.vPrice=raw_data.v_EV;
ddm_trace_data.vLotteryBuy=raw_data.v_EV;
ddm_trace_data.vLotterySell=raw_data.v_EV;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);
ddm_trace_data.zBuy=exp(raw_data.z_Intercept_trans)./(1+exp(raw_data.z_Intercept_trans));
ddm_trace_data.zDelta=raw_data.z_C_Role__T_2_;
ddm_trace_data.zSell=ddm_trace_data.zBuy+ddm_trace_data.zDelta;
ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

%% DDM4ae: without response tendency

raw_data=readtable([root_path,'\model\hddm\DDM4ae\DDM4ae_trace.csv']);

ddm_trace_data=table;
ddm_trace_data.t=raw_data.t;
ddm_trace_data.aBuy=raw_data.a_Intercept;
ddm_trace_data.aDelta=raw_data.a_C_Role__T_2_;
ddm_trace_data.aSell=raw_data.a_Intercept + raw_data.a_C_Role__T_2_;
ddm_trace_data.vInterceptBuy=raw_data.v_Intercept;
ddm_trace_data.vInterceptDelta=raw_data.v_C_Role__T_2_;
ddm_trace_data.vInterceptSell=raw_data.v_Intercept + raw_data.v_C_Role__T_2_;
ddm_trace_data.vPrice=raw_data.v_Price;
ddm_trace_data.vLotteryBuy=raw_data.v_Lottery_C_Role__1_;
ddm_trace_data.vLotterySell=raw_data.v_Lottery_C_Role__2_;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);
ddm_trace_data.zBuy(:)=0.5;
ddm_trace_data.zDelta(:)=0;
ddm_trace_data.zSell(:)=0.5;
ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

%% DDM4af: without a

raw_data=readtable([root_path,'\model\hddm\DDM4af\DDM4af_trace.csv']);

ddm_trace_data=table;
ddm_trace_data.t=raw_data.t;
ddm_trace_data.aBuy(:)=1;
ddm_trace_data.aDelta(:)=0;
ddm_trace_data.aSell(:)=1;
ddm_trace_data.vInterceptBuy=raw_data.v_Intercept;
ddm_trace_data.vInterceptDelta=raw_data.v_C_Role__T_2_;
ddm_trace_data.vInterceptSell=raw_data.v_Intercept + raw_data.v_C_Role__T_2_;
ddm_trace_data.vPrice=raw_data.v_Price;
ddm_trace_data.vLotteryBuy=raw_data.v_Lottery_C_Role__1_;
ddm_trace_data.vLotterySell=raw_data.v_Lottery_C_Role__2_;
ddm_trace_data.wLotteryBuy=ddm_trace_data.vLotteryBuy./(-ddm_trace_data.vPrice);
ddm_trace_data.wLotterySell=ddm_trace_data.vLotterySell./(-ddm_trace_data.vPrice);
ddm_trace_data.zBuy=exp(raw_data.z_Intercept_trans)./(1+exp(raw_data.z_Intercept_trans));
ddm_trace_data.zDelta=raw_data.z_C_Role__T_2_;
ddm_trace_data.zSell=ddm_trace_data.zBuy+ddm_trace_data.zDelta;
ddm_trace_data.val=ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy;
ddm_trace_data.res=ddm_trace_data.zSell - ddm_trace_data.zBuy;
ddm_trace_data.vInterceptBias=ddm_trace_data.vInterceptSell - ddm_trace_data.vInterceptBuy;
ddm_trace_data.aBias=ddm_trace_data.aSell - ddm_trace_data.aBuy;
ddm_trace_data.vallog=log(ddm_trace_data.wLotterySell ./ ddm_trace_data.wLotteryBuy);

% mean and 95% credible intervals
A=mean(table2array(ddm_trace_data));
B=prctile(table2array(ddm_trace_data),2.5);
C=prctile(table2array(ddm_trace_data),97.5);
rowNames = {'mean','2.5p','97.5p'};
colNames = {'t','aBuy','aDelta','aSell','vInterceptBuy','vInterceptDelta','vInterceptSell','vPrice','vLotteryBuy','vLotterySell','wLotteryBuy','wLotterySell','zBuy','zDelta','zSell','val','res','vInterceptBias','aBias','vallog'};
ddm.group.post.describe = array2table([A;B;C],'RowNames',rowNames,'VariableNames',colNames);
disp(ddm.group.post.describe);

%%










%% DIC
dic_raw=readtable([root_path,'\model\hddm\DDM4a\DDM4a_dic.csv']);
dic.DDM4a=dic_raw.dic(1);
dic_raw=readtable([root_path,'\model\hddm\DDM2a\DDM2a_dic.csv']);
dic.DDM2a=dic_raw.dic(1);
dic_raw=readtable([root_path,'\model\hddm\DDM3a\DDM3a_dic.csv']);
dic.DDM3a=dic_raw.dic(1);
dic_raw=readtable([root_path,'\model\hddm\DDM4b\DDM4b_dic.csv']);
dic.DDM4b=dic_raw.dic(1);
dic_raw=readtable([root_path,'\model\hddm\DDM4\DDM4_dic.csv']);
dic.DDM4=dic_raw.dic(1);
dic_raw=readtable([root_path,'\model\hddm\DDM0a\DDM0a_dic.csv']);
dic.DDM0a=dic_raw.dic(1);

%%










%% Surplus difference
sub_data_beh=readtable([root_path,'\data\subject\subject_behavior.csv']);
sub_data_beh.SurplusDiff=0-sub_data_beh.WTP-sub_data_beh.WTA;
[r,p]=corr(sub_data_beh.SurplusDiff,sub_data.aBias);
[r,p]=corr(sub_data_beh.SurplusDiff,sub_data.ValBias);
[r,p]=corr(sub_data_beh.SurplusDiff,sub_data.ResBias);
[r,p]=corr(sub_data_beh.SurplusDiff,sub_data.vInterceptBias);

%%
toc;
